Real-time dynamics in Quantum Impurity Systems: A Time-dependent Numerical 

Renormalization Group Approach 
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We develop a general approach to the nonequilibrium dynamics of quantum impurity systems for 
arbitrary coupling strength. The numerical renormalization group is used to generate a complete 
basis set necessary for the correct description of the time evolution. We benchmark our method 
with the exact analytical solution for the resonant-level model. As a first application, we investigate 
the equilibration of a quantum dot subject to a sudden change of the gate voltage and external 
magnetic field. Two distinct relaxation times are identified for the spin and charge dynamics. 
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Introduction. The investigation of real-time dynam- 
ics in quantum impurity systems (QIS) is of prime im- 
portance for our understanding of dissipation and deco- 
herence in qubits and electronic transport through nano- 
devices. Such systems consist of a mesoscopic subsystem 
or device, interacting with an infinitely large environment 
of noninteracting particles. Typical examples are quan- 
tum dots, qubits or biological donor-acceptor molecules. 
These are usually modeled either by a spin-Bose model, 
where a two-level system or qubit is interacting with 
a Bosonic bath 1], or by a set of atomic orbitals cou- 
pled to electronic leads, as in ultra-small quantum dots 
and single-electron transistors 0. Recently, Elzerman 
et al. reported the usage of gate-voltage pulses for a 
single-shot read out of the spin configuration of a single- 
electron transistor in a finite magnetic field. Such devices 
may be suitable for quantum computing, as they combine 
a long-lived quantum state with an easy read-out proto- 
col furnished by the coupling to the environment. 

In contrast to equilibrium conditions, the understand- 
ing of the real-time evolution of many-particle quantum 
systems is still at its infancy. Over the past 40 years, 
the Keldysh technique Q has proven to be the most suc- 
cessful approach to nonequilibrium (NEQ) dynamics, as 
it furbishes a perturbative expansion of the density op- 
erator. In general, however, perturbation theory fails in 
QIS due to the infra-red divergences caused by degenera- 
cies on the impurity. Recently, there has been significant 
progress in extending the density matrix renormaliza- 
tion group (DMRG) to time-dependent ID quantum sys- 
tems HIE GUI- The so-called TD-DMRG works well for 
finite-size systems and short times, but has an accumu- 
lated error proportional to the time elapsed. This makes 
the TD-DMRG unsuitable at present for tackling long 
time scales, most notably the exponentially long time as- 
sociated with the development of the Kondo effect 

In this Letter, we develop an alternative approach to 
the real-time dynamics of QIS based on Wilson's numeri- 
cal renormalization-group (NRG) method 0] • The NRG 
is a very powerful and accurate tool for calculating equi- 



librium properties of arbitrarily complex quantum impu- 
rities. We extend the approach to a certain class of time- 
dependent problems where a sudden perturbation, e.g., a 
gate- voltage pulse Q, is applied to the impurity at time 
t = 0. An earlier attempt [llj at a time-dependent NRG 
(TD-NRG) was faced with conceptual problems how to 
adequately couple the low- and high-energy scales. Here 
we resolve these problems by combining the NRG with 
Feynman's concept of a reduced density matrix [l^ . Ilflj , 
and by implementing a suitable resummation procedure 
for tracking all states that contribute to the time evo- 
lution of the system. We establish the accuracy of our 
approach by comparison to the exact analytical solution 
of the resonant-level model, and apply it to the single im- 
purity Anderson model (SIAM). Two distinct time scales 
are found for the spin and charge dynamics in the SIAM. 
The process of spin relaxation is shown to be sensitive to 
the initial conditions imposed on the system. 

Formulation. As our approach applies to arbitrary 
QIS, we first formulate it in general terms before turn- 
ing to concrete examples. The Hamiltonian of QIS is 
generally given by W = TL ba th + H 



imp 



Tin 



where 



Tibath models the continuous bath, TLimp represents the 
decoupled impurity, and Ti m i X describes the coupling be- 
tween the two subsystems. The entire system is assumed 
to be in thermal equilibrium for t < 0, at which point 
a sudden perturbation A7i is switched on: 7i(t > 0) = 
H l + AH = HJ . For t > 0, the density-matrix operator 
evolves according to pit > 0) = e~ ltut p eq e ltut , where 
p eq = e-^/Trie-W}. 

The key ingredient in the NRG is a logarithmic dis- 
cretization of the continuous bath, controlled by the pa- 
rameter A > 1 10]. The Hamiltonian is mapped onto 
a semi-infinite chain, where the Ath link represents an 
exponentially decreasing energy scale Dm ~ A~ N ^ 2 [Icj . 
Using this hierarchy of scales the sequence of finite-size 
Hamiltonians TLn is solved iteratively, truncating the 
high-energy states at the conclusion of each step to main- 
tain a manageable number of states. The reduced basis 
set oi Hn so obtained is expected to faithfully describe 
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FIG. 1: The full NRG chain of length N is divided into a 
sub-chain of length m and the "environment" R m ,N- 



the spectrum of the full Hamiltonian on a scale of Dn, 
corresponding to the temperature T/v ~ E>n 10]. 

Out of equilibrium, one can no longer settle with the 
reduced basis set of at temperature Tjv, as the NEQ 
dynamics involves all energy scales exceeding Tjy. A com- 
plete basis set of the Fock space Tn of Hn is required. 
One possible choice of a complete basis set is given by 
Tn — {\l,e;m) tr }, where I labels the NRG eigenstates 
of Ti.^ discarded at iteration to = 1, • • • , N, and e is a 
tuple of quantum numbers labeling the "environment" 
degrees of freedom on the remaining chain tail R m ,N, 
see Fig. For to — N, all NRG eigenstates are consid- 
ered discarded. Restricting attention to local operators 
O which do not act on the "environment" degrees of free- 
dom, the time evolution of (O) is formally given by 



(o)(t) 



E 



^ tr(l,e;m\d\l',e';m')tr 



tril^e'im'le-^p^e^lheim}^ . (1) 



The sum over to, ml in Eq. (yi divides into three parts: 
m > ml , to = to', and m < ml . Consider m > m! . 
Introducing the states {\l, e; m)k P }, where I now labels 
the retained NRG eigenstates of H{ n and e records the 
"environment" degrees of freedom, one has the identity 



^ \l',e';m') tr tr(l',e';m'\ 



m'.l'.e' 



\l, e; m)k P k P (l, e; m\. 



(2) 



A similar identity applies to ml > m upon interchanging 
the roles of m and ml . Inserting Eq. (J2J) and its to <-> ml 
counterpart into Eq. Q one obtains 



(0)(t) = 



N trun 

EE 

m— r,s 



Equation © is one of the central results of this paper. 
In contrast to the TD-DMRG d H %M and a previous 
attempt at a time-dependent NRG TjJ , all states of the 
finite-size Fock space are retained, and all energy scales 
are explicitly taken into account. No basis set reduction is 
imposed. To mimic the relaxation in an infinite-size sys- 
tem we (i) average Eq. © over N z different realizations 
of the NRG discretization using Oliveira's z-trick 
and (ii) introduce a scale-dependent damping factor, 



i(E™-E™)t 



e i{E™-E T )t e - 



xD m t 



(5) 



for all E™ ^ E™ in Eq. J5J| . Here a is a constant of 
order unity, representing a Lorentzian broadening of the 
NRG levels. Below we present results with and without 
the additional damping factor a. 

Algorithm. To implement the TD-NRG at temper- 
ature T, one first selects a chain length N such that 
T ks Tjy. Two simultaneous NRG runs are then per- 
formed, one for HI and another for 7i* , to obtain (i) 
the equilibrium density matrix p eq , (ii) the NRG eigen- 
energies and eigenstates of W and H l , and (iii) the 
overlap matrices S r y (to) between the eigenstates of Til 
and TL^ at each iteration m. Using a recursion relation 
p red (m) p red (m — 1), the reduced density matrix of 
Eq. Q is evaluated with respect to the NRG eigenstates 
ofW, and then rotated to the NRG eigenstates of TL^ us- 
ing the overlap matrices S r>r i(m). The expectation value 
(O) (t) follows from evaluation of Eq. © . 

Benchmark. As a critical test of our method, we first 
apply it to the resonant-level model (RLM), describing 
the hybridization of a localized level d) with a band of 
spinless conduction electrons c\: 

H = E £kC l Ck + E d^d + Vj2{^ c k + c{d}. (6) 



Choosing fid = d^d as the observable O, we consider 
a stepwise change in the energy of the level: Ed(t) = 
6(-t)E^+6(t)E d . In the wide-band limit, n d (t) = (n d )(t) 
can be solved exactly in closed analytical form using the 
Keldysh formalism. For T = 0, the solution features 
an exponential decay from the initial equilibrium occu- 
pancy of to the new equilibrium occupancy of W 
with two decay rates T and 2T. Terms decaying at rate 



where O™ = (r, e; to|0|s, e; to) is independent of e, and £ show an additional modulation of frequency \Ej \ 1 
p r s ed is the reduced density matrix |l2l Il3| 



red 
rs.r 



(to) 



^2(s,e;m\pe q \r,e;m}. 



(4) 



Here we have used the conventional NRG approximation 
H f N \l,e;m) w £ ; m |Z,e;m), where Ef 1 is the NRG eigen- 
energy of corresponding to the eigenstate I . The 
restricted sum X)r "™ m ^1- El implies that at least one 
of the states r and s is discarded at iteration to. 



Here T = irp^V 2 is the hybridization width, po being the 
conduction-electron density of states at the Fermi energy. 

Figure |21 shows a comparison of the TD-NRG with the 
exact analytical solution for rid{t), in response to sudden 
change in the energy of the level from E% = to E\ ^ 0. 
Here we have used a symmetric box density of states 
p{e) = p 6(D - |e|) for the TD-NRG El and averaged 
over N z = 32 equally spaced values of z [14j . No damping 
factor was introduced. To facilitate comparison with the 
wide-band limit, a large bandwidth D/T — 500 was used. 
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FIG. 2: Time-dependent occupancy of the RLM for (a) E® = 
0, T -> and three different values of E^/T = -1, -2, -3, 
and (b) E° d = 0, E\ = —F and four different temperatures 
T/T = 0.1, 0.5, 1, 5. Solid lines are the results of the TD-NRG 
with A = 1.3, D/T = 500, N s = 1000, N z = 32, and a = 
(no damping factor). Dashed lines are the exact solution in 
the wide-band limit. Right arrows mark the saturated NRG 
occupancies for a > 0. Inset to (a): Progressive calculation 
of nd(0) as a function of the number of backward iterations 
N-m, for Eh = -2T. 



As seen in Fig. Ota), the TD-NRG accurately repro- 
duces the exact T = solution on all time scales up to 
moderately large E\. For E\ = —T, the maximal error 
is less than 1%. For E\ = -2T it is less than 2.5%. The 
short-time behavior is well described by the TD-NRG for 
all E\. In particular, n c i{t — > + ) coincides with the ini- 
tial equilibrium NRG occupancy of TC independent of E° 
and Ej. As demonstrated in the inset to Fig.|2a), sum- 
mation over the states of all NRG iterations is essential 
for recovering nd(0). One cannot settle with just a single 
NRG iteration as in the equilibrium case. 

Remarkably, the TD-NRG retains the same accuracy 
in Fig.|2Ia) up to arbitrarily long times. Deviations from 
the new equilibrium occupancy of Hj are only slight at 
large t, and are further reduced by decreasing A, due to 
elongation of the effective chain length corresponding to 
the energy scale E\. The long-time deviations in rid(t) 
stem from a difference between the time-evolved density 
matrix in a finite-size system and the equilibrium density 
matrix p oc e~^ w , approached only for A — > 1 + . 

The TD-NRG can also be applied to finite tempera- 
tures, see Fig.|2Jb). For moderate E\, there is excellent 
agreement with the exact curves up to t several times 
larger than 1/T, at which point finite-size oscillations 
develop. The latter oscillations are greatly reduced by 
averaging over the different z's, to the extent that they 
appear damped at longer times. Note, however, that all 
curves in Fig.[3fb) have practically saturated at their new 
equilibrium values prior to the appearance of oscillations. 

We stress that no explicit damping factor was utilized 



in Fig. [21 A nonzero a depresses the oscillations in n<i(t), 
which saturates at the new steady-state values indicated 
by arrows in Fig. 0a) . The difference 5nd between the 
new equilibrium occupancy of HJ and the steady-state 
value reached serves as a sharp criterion for the accuracy 
of the TD-NRG. The smaller Sn d is the better the TD- 
NRG performs on all time scales. This provides one with 
a valuable estimate for the accuracy of the TD-NRG in 
cases where no other reference exists. 

Single impurity Anderson model. The first nontrivial 
model which lacks an exact solution out of equilibrium is 
the single impurity Anderson model (SIAM), describing 
a spinful resonant level with an on-site repulsion U: 

n = $>** c L < w+ E [E d (t)-^H(t)]4d a 

V(t)J2{4ck, + cld a } (7) 



+Uh d n d + 



k.a 



with n d = S a d a . This Hamiltonian is commonly used 
to model a single Coulomb-blockade resonance in ultra- 
small quantum dots. To probe the relaxation both in the 
spin and charge sectors, we consider in the following a 
simultaneous stepwise change in the energy of the level, 
Edit) = 9{-t)E° d + 9{t)E\, and in the local magnetic 
field, H(t) = 9(-t)H + 6(t)Hi. We also permit a sud- 
den change in the hybridization width T(t) = irpoV 2 (t): 
r(t) = 0{— t)To + 8(t)Ti. As our prime interest is in the 
relaxation to a new Kondo state, we focus hereafter on 
Hi = and only briefly comment on nonzero Hi . 

In the spirit of the experiment of Elzerman et al. 0], 
we begin with a spin-polarized impurity, where the lower 
spin state is degenerate with the empty-orbital state: 
E d - #0/2 = 0. At t = 0, both the local magnetic field 
is switched off and the energy of the level is shifted to 
E\ = —{7/2. The resulting time evolutions of the impu- 



rity occupancy nd(t) 



f)(t) and magnetization 



if) (t) are plotted in Fig. for T -> 0, 



Ho = Ti, and different values of U. Two separate initial 
conditions are considered: Fo = Ti and To = 0, repre- 
senting different relaxations to the strong-coupling fixed 
point. The former corresponds to relaxation from a spin- 
polarized mixed-valent state, while the latter describes 
evolution from a free- impurity fixed point. 

Two distinct time scales are clearly visible in Fig.|3 for 
the relaxation of spin and charge. Similar to the RLM, 
n>d{t) equilibrates on a time scale £ c h cx 1/Fi, and de- 
velops Rabi-type oscillations for \E\\ > Ti. For t 3> £ c h, 
the occupancy must saturate at its new equilibrium value 
rid = 1< This behavior is well captured by the TD-NRG, 
attesting to its accuracy. The long-time deviations Snd 
are less than 0.015 for Tq = Ti and less than 0.04 for 
Tq = 0. Importantly, Srid is notably smaller for the SIAM 
than for the RLM with comparable E\. 

Contrary to the dynamics of nd(t), the spin magneti- 
zation S z (t) equilibrates on a much longer time scale t sp . 
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FIG. 3: Time- dependent occupancy rid(t), Figs, (a) and (b), 
and magnetization S z (t), Figs, (c) and (d), of the SIAM with 
T -> 0, D/Ti = 20, Ho = 2E° d = Ti, Hi = 0, E\ = -(7/2, 
and six different values of U/Yi = 2, 4, 6, 8, 10, 12. Here To = 
Ti in Figs, (a) and (c), whereas Fo = in Figs, (b) and (d). 
NRG parameters: A = 1.8, N a = 1000, N z = 16, and a = 0.1. 



To elucidate the strong U dependence of t sp , we replot- 
ted the magnetization curves of Fig. [3] versus t/tK with 
t K = 1/T K , see Fig. 0] Here T K is the Kondo temper- 
ature of H f , defined by T KX {T K ) = 0.07 \vj^ ix being 
the impurity susceptibility of Ti^). While the initial dy- 
namics of S z (t) is governed by i c h oc 1/ri, the eventual 
relaxation is clearly governed by the Kondo time scale: 

t sp oc 1/T K H3- This is best seen for r o = 0, Fig.Efb), 
where all curves gradually collapse onto a single master 
curve, provided a stable local moment was formed on the 
level. For Y$ = T±, there is no universality in t/tK- In 
fact, S[t) differs in form from that of r = 0, revealing 
a sensitivity of spin relaxation to the initial conditions 
imposed on the system. A nonzero H\ suppresses the 
Kondo effect. Accordingly, the spin relaxation time t sp 
is reduced, approaching i c h for Hi ^> Tk (not shown). 

Summary and discussion. A powerful new approach 
to the noncquilibrium dynamics of QIS was presented 
based upon Wilson's NRG. All states of the Wilson chain 
are summed over in this scheme, and no basis set reduc- 
tion is imposed. We established the accuracy of the ap- 
proach on all time scales by comparison to the exact so- 
lution of the RLM, and demonstrated its ability to access 
exceedingly long time scales by investigating the equili- 




FIG. 4: The impurity magnetization curves of Fig. |3] replot- 
ted versus t/tK with tK = 1/Tk- Here To = Ti in Fig. (a), 
while Tq = in Fig. (b). 



bration of the SIAM. The spin relaxation time t sp of the 
SIAM was shown to be governed by the Kondo time scale 
tx = 1/Tk, for T -C Tk- However, there is no univer- 
sality in t/tK- Spin relaxation is sensitive to the initial 
conditions imposed on the system. Although our specific 
examples were focused on Fermionic baths, the TD-NRG 
can equally be applied to Bosonic baths [18j . This should 
allow for an accurate investigation of the time evolution 
of a two-level system [l| , as relevant to quantum comput- 
ers, and of biological donor-acceptor molecules. Finally, 
we expect that our approach can be generalized to contin- 
uous time dependences (such as ac drives) by discretizing 
the time axis, and evaluating Eq. J2J consecutively using 
the discrete set of Hamiltonians 7i*{tj). This will pro- 
vide a very general approach to the real-time dynamics 
of quantum-impurity systems at any temperature T. 
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der project no. HHB000. A.S. was supported in part by 
the Centers of Excellence Program of the Israel Science 
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